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ABSTRACT 

On 1998 November 14, Saturn and its rings occulted the star GSC 0622-00345. The occultation latitude 
was 55.°5 S. This paper analyzes the 2.3 fim light curve derived by Harrington & French. A fixed-baseline 
isothermal fit to the light curve has a temperature of 140 ± 3 K, assuming a mean molecular mass of 2.35 
AMU. The thermal profile obtained by numerical inversion is valid between 1 and 60 //bar. The vertical 
temperature gradient is >0.2 K km"' more stable than the adiabatic lapse rate, but it still shows the alternating- 
rounded-spiked features seen in many temperature gradient profiles from other atmospheric occultations and 
usually attributed to breaking gravity (buoyancy) waves. We conduct a wavelet analysis of the thermal profile, 
and show that, even with our low level of noise, scintillation due to turbulence in Earth's atmosphere can 
produce large temperature swings in light-curve inversions. Spurious periodic features in the "reliable" region 
of a wavelet amplitude spectrum can exceed 0.3 K in our data. We also show that gravity-wave model fits to 
noisy isothermal light curves can lead to convincing wave "detections". We provide new significance tests for 
localized wavelet amplitudes, wave model fits, and global power spectra of inverted occultation light curves by 
assessing the effects of pre- and post-occultation noise on these parameters. Based on these tests, we detect 
several significant ridges and isolated peaks in wavelet amplitude, to which we fit a gravity wave model. We 
also strongly detect the global power spectrum of thermal fluctuations in Saturn's atmosphere, which resembles 
the "universal" (modified Desaubies) curve associated with saturated spectra of propagating gravity waves on 
Earth and Jupiter. 

Subject headings: atmospheric effects — methods: statistical — occultations — planets and satelUtes: atmo- 
spheres — planets and satellites: individual (Saturn) — waves 



1. INTRODUCTION 

Earth-based occultations remain an attractive method for 
measuring the thermal profile in the l-100-/xbar region of a 
planetary atmosphere. Many such profiles for Saturn were 
recorded during the 28 Sgr occultation of 1989 July 3, which 
sampled the equatorial region from 6.°6 N-15.°2 S latitude 
(Hubbard et al. 1997). There is a single profile for the north 
polar region (Cooray et al. 1998, 82.°5-85° N), and a north- 
ern low-latitude profile from the same event (French et al. 
1999, 19.° 16 N). Saturn's central flash probes much deeper, 
around 2.5 mbar; Nicholson et al. (1995) obtained IR images 
of the flash during the 28 Sgr event, from which they inferred 
the zonal wind profile of the sampled latitudes along Saturn's 
limb. 

Occultations observed by a spacecraft near a giant planet 
probe the troposphere from the cloud deck (^1 bar) to the 
mbar level at radio and infrared wavelengths. They probe the 
upper stratosphere and thermosphere (<1 //bar) in the ultravi- 
olet. Earth-based visual and infrared occultations measure the 
thermal structure of the intervening mesosphere and strato- 
sphere regions, which are not well sampled by spacecraft ex- 
periments. For Saturn, the Pioneer radio (Kliore et al. 1980; 
Lindal et al. 1985) and Voyager 2 extreme ultraviolet solar 
and stellar observations (Smith et al. 1983) sensed the equa- 
torial region only. The Voyager 1 radio occultation sensed 75° 

"^ Visiting Astronomer at the Infrared Telescope Facility, which is oper- 
ated by the University of Hawaii under Cooperative Agreement no. NCC 
5-538 with the National Aeronautics and Space Administration, Science 
Mission Directorate, Planetary Astronomy Program. 



S (Tyler et al. 1981), while the Voyager 2 radio occultations 
sensed 36.°5 N and 31° S (Tyler et al. 1982). The Cassini ra- 
dio experiment has performed a number of radio occultations 
in the equatorial region as well as at middle and high latitudes 
of Saturn (Nagy et al. 2006; Kliore et al. 2009). The Cassini 
Ultraviolet Imaging Spectrograph stellar occultation probed 
the upper atmosphere at 40° S and 66° N (Shemansky 2008), 
and the Cassini Composite Infrared Spectrometer mapped 
Saturn's thermal atmospheric emission, resulting in temper- 
ature maps for both hemispheres for pressures ranging from 
0.1 mbar to about 700 mbar (Flasar et al. 2005). The Cassini 
Visual and Infrared Mapping Spectrometer also has the ca- 
pability of observing spectrally resolved near-infrared stellar 
occultations by Saturn's atmosphere (Brown et al. 2004). 

Temperature profiles for a variety of atmospheres from 
both occultations and in situ observations show quasi-periodic 
structures that are usually attributed to propagating waves. 
Waves have been reported on Venus (Hinson & Jenkins 1995), 
Earth (Fritts & Alexander 2003, and references therein). Mars 
(Creasey et al. 2006; Fritts et al. 2006), Jupiter (French & 
Gierasch 1974; Young et al. 1997; Raynaud et al. 2003, 2004; 
Young et al. 2005), Saturn (Cooray et al. 1998; Fouchet et al. 
2008), Titan (Sicardy et al. 1999), Uranus (Young et al. 2001), 
Neptune (Roques et al. 1994), and Pluto (Person et al. 2008; 
Hubbai-d et al. 2009; Toigo et al. 2010). Both the behavior 
of individual waves and the form of wave power spectra can 
reveal properties of the underlying atmosphere. For example, 
the forcing, propagation, and dissipation of the waves both 
contribute to and depend on the sources and sinks of energy 
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Table 1 

Isothermal Fit Results 



Parameter 



Baselines Free 



Baselines Fixed 



Half-light time (UTC) 

Full flux 

Background 

Scale height, H (km) 

Temperature, T (K) 



11:18:47.14 ±0.18 

1.0157 ±0.0028 

-0.0150 ± 0.0013 

50.7 ± 1.1 

159 ±4 



11:18:46.59 ±0.17 

1 (not fit) 

(not fit) 

44.7 ± 0.9 

140 ± 3 



in the atmosphere, the background thermal state, and eddy 
and molecular diffusion. 

On 1998 November 14, Saturn and its rings occulted GSC 
0622-00345, as predicted by Bosh & McDonald (1992). We 
obtained a light curve for atmospheric immersion, based on 
infrared imaging observations at the NASA Infrared Tele- 
scope Facility (IRTF) on Mauna Kea, HI. The high signal- 
to-noise ratio (S/N) allowed us to determine the vertical tem- 
perature profile of Saturn's stratosphere at 55. °5 S latitude, a 
region not sampled by previous stellar occupation observa- 
tions (see Figure 1 and Table 1 of Harrington & French 2010, 
hereafter Paper I). 

Paper I presents the light curve and describes the new meth- 
ods used to acquire and derive it. This paper presents the sci- 
entific analysis of the light curve. Subsequent sections cover 
isothermal model fits, numerical inversions to derive the ther- 
mal profile, noise tests, local and global wavelet spectrum 
analysis, a gravity-wave model based on wavelet reconstruc- 
tion, exploration of the "universal" power spectrum of gravity 
waves, discussion of the global power spectrum, and our con- 
clusions. For each analysis, we present new significance tests 
that determine the effects of real (non-Gaussian) noise. 

2. ISOTHERMAL FITS 

Table 1 presents the results of isothermal model fits to the 
light curve with free and fixed baselines. The free baselines 
flank the calculated values by over 1.5% of full flux, a huge 
deviation given the high accuracy of the baseline determina- 
tion. The second fit fixes the baselines at their calculated val- 
ues. For the latter case, the derived scale height is very close 
to that in the good region of the inversion presented below. 

Isothermal fits can give approximate light-curve parameters 
that do not depend on many assumptions and that are unique. 
However, the vastly different values for temperature, T, and 
scale height, H, between the two cases and the poorly fit base- 
lines in the free-baselines case indicate that isothermal models 
do not approximate this atmosphere well. Both models have 
long, nonzero tails that stand well above the data (see Figure 
2 of Paper I). The light curve also contains many spikes with 
amplitudes that are many times the noise level. The spikes 
and the low-valued tail are features of the observations that 
are not in the model but that strongly influence where the fit 
falls. Introducing a temperature gradient to the model might 
improve the fit, but that does not address the numerous spikes. 
Inversion is thus the proper analytic approach. 

Both here and in the inversion that follows, we use the mean 
molecular mass to convert H toT. The He/H2 volume mixing 
ratio for Saturn determined by Conrath & Gautier (2000) is 
in the range 0.11-0.16. We adopt a value of 0.135 and a CH4 
mixing ratio of (4.7 ±0.2) x 10"^ (Fletcher et al. 2009), result- 
ing in a mean molecular mass of 2.35 AMU. The uncertainty 
in He/H2 dwarfs the uncertainty in any other constituent. The 
value of 2. 1 35 AMU used by Hubbard et al. ( 1 997) and many 
prior workers is much smaller. One must be careful to ad- 



just temperatures and adiabatic lapse rates to the same mean 
molecular mass when making comparisons. 

3. TEMPERATURE PROFILE 

We applied an Abel transform to the normalized light curve 
under the usual assumptions that the atmosphere is radially 
symmetric and that ray crossing is not substantial (French 
et al. 1978; in grazing occultations, unlike ours, these assump- 
tions may be violated, see Cooray & Elliot 2003). This pro- 
duced T, pressure (P), and number density (n) as a function 
of height above the half-light level (z). Table 2 and Figure 1 
show our parameters and results. We rebinned the light curve 
prior to the inversion into time intervals corresponding to ver- 
tical atmospheric layers of equal thickness (1 km). This is 
high enough resolution to preserve the intensity spikes in the 
light curve and is roughly comparable to the 1.2 km Fresnel 
scale (Paper I). Saturn is very oblate, so at any point on the 
surface the radius of curvature depends on both the latitude 
and the direction being considered (e.g., north-south versus 
east-west). Using the radius of curvature along the line of 
sight at the half-light latitude accounts for the planet's oblate- 
ness. The derived vertical refractivity profile is proportional 
to the density profile. The T and P profiles then result from 
integrating the hydrostatic equation and applying the ideal gas 
law (French et al. 1978). 

The inversion process determines the unique vertical tem- 
perature profile such that a forward model of the occultation 
would reproduce the observed light curve exactly, subject to 
assumptions just mentioned. Unfortunately, the inversion can- 
not distinguish between noise and atmosphere-induced varia- 
tions in stellar intensity. Both random and systematic noise in 
the observations can thus seriously affect the derived tempera- 
ture profile (see Elliot et al. 2003 for an extensive review, and 
references cited therein). Initially, the photometric noise in 
the upper baseline completely overwhelms the actual refrac- 
tive effects of the tenuous upper atmosphere. The noise intro- 
duces ^ 100 K swings in the upper part of the derived thermal 
profile, and these unphysical swings bias the results at the on- 
set of the inversion. In deeper layers, the uncertainties associ- 
ated with this initial condition are less significant. Eventually, 
refractive defocussing of the starlight is so large that the noise 
in the lower baseline of the light curve dominates the faint 
signal. This results in unphysical temperature variations and 
trends at the deepest atmospheric layers probed during the oc- 
cultation (see, e.g., Raynaud et al. 2004). As a rough guide, 
French et al. (1978) showed that, for high-quality. Earth-based 
stellar occultations, the valid region of the derived tempera- 
ture profile ranges from about 0.5 to -3.5 H above and below 
the half-light level. 

Several strategies have been adopted to minimize the effects 
of the unstable initial condition of the inversion process. El- 
liot & Young (1992) and Elliot et al. (2003) explored a range 
of models that fit the upper part of the light curve assuming a 
power-law dependence of temperature with radius. This is an 
extension of the strategy developed by French et al. (1978) of 
replacing the upper part of the light curve by the best isother- 
mal fit to that restricted part of the data. In the absence of de- 
tailed knowledge of the upper stratospheric temperature struc- 
ture above the occultation region, we adopt the minimalist as- 
sumption that the atmospheric region sounded by the upper 
part of the light curve is isothermal. Effectively, we assume 
that the upper atmosphere above the inversion's reliable re- 
gion does not have large-scale temperature fluctuations, and is 
comparable in mean temperature to the valid region. A com- 
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Figure 1. Family of thermal profiles derived by numerical inversion. The light curve inverted for each curve has an isothermal cap that ends at the flux level 
indicated in the key. We used the 75% curve (green line in the middle of the group) for subsequent analyses. The data are plotted vs. pressure (left axis), with the 
approximate corresponding altitude scale given on the right axis. Left: temperature profiles. The dashed line shows the temperature of our isothermal fit. Middle: 
the same data after removal of a linear fit to each profile. This shows that the shape and amplitude of the small-scale temperature fluctuations is consistent among 
the profiles, that there are veiy small differences in their altitudes, and that the deviations from a linear thermal profile are no more than about ±4 K. Right: 
vertical temperature gradient vs. altitude, derived from the left panel's data. The atmosphere is statically stable over this altitude range, since the temperature 
gradient is separated from the adiabatic lapse rate (dashed line). Nonetheless, we still see alternating rounded and spiked structures, which are seen in many other 
profiles where the vertical temperature gradient approaches adiabatic. 



parison of 28 Sgr and Voyager UV stellar occultation obser- 
vations bears this out for Saturn (Hubbard et al. 1997). The 
stability of the hybrid-light-curve inversion depends on the 
length of this isothermal "cap." If the cap extends only from 
the upper baseline to the 99% level (in units of normalized 
stellar flux), then the noise in the subsequent upper part of the 
observed light curve will still produce spurious temperature 
variations at the onset of the inversion. As the cap length in- 
creases, the inversion stabilizes, eventually contaminating the 
inversion's valid region. The optimal cap is large enough to 
give a stable inversion but ends above the valid region. 

The left panel of Figure 1 shows the thermal structure de- 
rived from a suite of light-curve inversions with seven differ- 
ent isothermal caps ending at 97.5% - 70% of the full stellar 
signal (we computed additional caps outside this range). The 
order of the curves reflects the sensitivity of the inversion to 
noise and atmospheric structure that immediately follow the 
isothermal cap. The degree of uncertainty in the inversion is 
shown by the spread between the profiles. We have used the 
75% cap in all subsequent analyses. Elliot et al. (2003) show 
that a 50% cap results in about a 3% error in the derived tem- 
perature at the onset of the underlying inversion region, for an 
isothermal light curve with S/N per scale height of 200 in the 
presence of white noise. Since our goal is to investigate ver- 
tical variations in temperature as well as to estimate the mean 
temperature, we have adopted a compromise of replacing only 
the upper 25% of the light curve by an isothermal model. 

The temperature variations are more clearly seen in the mid- 
dle panel of Figure 1, which shows the deviations of each of 
the profiles from a linear fit to that profile. The right panel 
shows the vertical temperature gradient for each inversion. 
The vertical dashed line corresponds to the adiabatic lapse 
rate. The atmosphere is locally stable against convection from 
1-60 fihai. Some occultation temperature gradient profiles 
(e.g.. Figure 10 of Raynaud et al. 2003, Figure 7 of Raynaud 
et al. 2004, and others cited therein) have oscillations with a 



Table 2 

Inversion Parameters and Results 



Description Value Units 


Comment 


Gravity at half-light, g H .06 m s"^ 
Refractivity at STP" 1.24x 10"^ 
Mean molecular mass 2.35 AMU 
Local rad. of curv. 64,307.9 km 
T/H 3.126 K km"' 
Adiabatic lapse rate, r -0.992 K km"' 


Conrath & Gautier (2000) 
Along line of sight 



" standard temperature and pressure 

rounded shape on the low AT /Az side and a narrow, spiked 
shape on the high side, and so does ours, particularly at depth. 
The asymmetric rounding has been attributed to gravity wave 
breaking as the profile's gradient approaches the adiabatic 
lapse rate. Saturn's adiabatic lapse rate is separated from 
the negative-side extrema in our profile's gradient by 0.2 K 
km"' everywhere, and generally by much larger amounts. One 
might expect that, as waves propagate vertically and increase 
in amplitude by virtue of energy conservation, they will even- 
tually become superadiabatic and result in wave breaking, but 
this is not always reflected in the retrieved vertical tempera- 
ture profiles. Similar results to ours have been found for Pluto 
(see Figure 4 of Young et al. 2008). The occultation profile in- 
tegrates over a large atmospheric path length, so there may be 
local gradient instabilities not seen in occultation inversions. 

4 NOISE TESTS 

It is not a simple matter to quantify in detail the effect of 
light-curve noise on the numerical inversions. The inferred 
structure at a given atmospheric level is contaminated by er- 
rors in the derived refractivity of all overlying levels. French 
et al. (1978) showed that the correlation length scale of the in- 
version process extends to well over a scale height above any 
given pressure level, as dictated by the width of the kernel in 
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Table 3 

Inversion Set Temperature Ranges and Power-law Fits to Global Wavelet Spectra 



Sample 


Temp. 


Temp. Amplitude 


Temp. Amplitude 


Mean-Norm T Pow 


Mean-Norm T Pow 


BG-Norm Temp. Pow 


BG-Norm Temp. Pow 


Set 


Range 


Mean 


Max 


Mean 


Max 


Mean 


Max 




Mean 


Const." Exp. 


Const. Exp. 


Const. Exp. 


Const. Exp. 


Const. Exp. 


Const. Exp. 




(K) 


(K) 


(K) 


(K) 


(K) 


(K) 


(K) 



9.85 
5.52 
1.60 
2.72 
0.84 



1.23x10- 



2.05 X 10-2 
3.71x10-^ 
1.79x10"^ 
3.27x10"^ 



-1.05 
-1.26 
-1.26 
-1.02 
-1.09 



4.12x10-2 
1.51x10"^ 
3.93x10"^ 
7.58x10'^ 



-1.23 
-1.22 
-0.99 
-1.07 



5.78 xUr*" 
2.76X Ur^ 
2.58x10"*^ 
1.52x10"^ 

7.19X ur^ 



-3.02 
-3.23 
-2.70 
-2.91 
-2.90 



1.16x10'*' 
4.96x10'^ 
6.16x10-^ 
3.32x10-^ 



-3.00 
-2.41 
-2.83 
-2.81 



5.44x10-' 
2.86x10-' 
2.69x10-* 
1.56x10"' 
7.40x10"- 



-3.05 

-3.22 1.25x10-*' 

-2.69 5.54x10-^ 

-2.91 6.33x10-^ 



-2.90 



3.42x10 



-8 



-2.96 
-2.40 
-2.83 
-2.81 



" Const, and Exp. refer to the constant and exponent (power) of the power-law fits (Section 7). 
'' Range is average range for different caps. Power law is for 75% cap. 



the integral equation for the Abel transform. Thus, even un- 
correlated white noise in the light curve results in correlated 
errors in the derived thermal profile. 

The problem is even more complex if the input light-curve 
noise is correlated, and almost any baseline drift on time 
scales longer than the exposure time will introduce signifi- 
cant correlation. Causes of such drifts include turbulence and 
waves in Earth's atmosphere ("seeing"), atmospheric trans- 
parency variations, and pointing drifts. The latter can be prob- 
lematic in spacecraft instruments with few spatial channels 
(e.g., few pixels) or with pointing-dependent sensitivity. 

Although previous investigators have explored the conse- 
quences of white noise to the mean temperature determined 
by inversion (e.g., French et al. 1978; Elliot et al. 2003), there 
has been no systematic study of the effects of noise on wave 
analyses, and little consideration of correlated noise. To as- 
sess the significance level of features in our derived profiles 
and subsequent analyses, we added both Gaussian and real 
noise to an isothermal model light curve. We created 25 re- 
alizations of Gaussian noise with the same standard deviation 
as our upper baseline. To eliminate the baseline uncertainty 
issue, we added this noise only below the 75% light level. 
For the real noise tests, we took a section of our upper base- 
line, removed a low-order polynomial, and repeated the sec- 
tion several times. The polynomial ensures that the sections 
have zero mean and meet without a discontinuity. We shifted 
the resulting data vector by three different amounts so that 
specific noise spikes would appear at three different locations 
along the synthetic light curve. We added this to the isother- 
mal curve starting at each of the cap levels to create 21 sets 
(seven caps times three shifts). 

Sources of scintillation noise include both the star and the 
residual from template subtraction (Paper I). To account for 
the decrease in the stellar contribution to scintillation, we cre- 
ated two additional sample sets, identical to those above but 
with the noise scaled by the normalized intensity of the noise- 
free isothermal model light curve. Most ground-based occul- 
tation data sets have similar noise levels on both baselines, in- 
dicating that the dominant noise source is residual planetary 
light and that the unsealed noise analysis is most appropriate 
for that case. However, reductions in the residual planetary 
signal are possible in the future, in which case the scaled noise 
analysis would apply. We label the four sets of noisy isother- 
mal profiles "RN" (real, normal), "RS" (real, scaled), "GN" 
(Gaussian, normal), and "GS" (Gaussian, scaled). 

We inverted all of these light curves and calculated the tem- 
perature ranges in the valid regions of each resulting profile. 
The first columns of Table 3 present the average of all the tem- 



perature ranges in each sample set and in the observations. 
Both types of noise induced large oscillatory structures into 
profiles that, without noise, should have been straight vertical 
lines, but real noise had a dramatically larger effect than Gaus- 
sian noise. Temperature ranges were typically 2 - 3 K for the 
GN set, although a single outlier had a 6.4 K range. Within 
a given shift of the RN set, the profiles and their ranges were 
mostly similar Ranges for the three shifts averaged 7.4, 5, 
and 4.5 K, indicating that the placement of individual light 
curve spikes strongly affected the results. The scaled sets had 
much smaller ranges, but otherwise behaved similarly to their 
unsealed brethren. 

The large-amplitude oscillations in these supposedly 
isothermal inversions generally have wavelengths of at least 
a scale height. Real noise is substantially worse than Gaus- 
sian noise of the same standard deviation, likely because of 
its red power spectrum. We assess the power spectrum of the 
noise and related error and significance issues in Section 7. 

5. WAVELET ANALYSIS 

The observed temperature profiles all show fluctuations 
with amplitudes as large as 4 K. Small-scale, quasi-periodic 
structures in atmospheric profiles are often interpreted as 
inertia-gravity waves (French & Gierasch 1974; Young et al. 
1997; Cooray et al. 1998; Raynaud et al. 2003, 2004). To 
investigate the wave nature of these structures in more de- 
tail, we computed the wavelet transform of the valid region. 
A wavelet transform gives the amplitude or power spectrum 
as a function of atmospheric depth. Torrence & Compo 
(1998) provide a quantitative and accessible wavelet tutorial 
with software; see Acknowledgements for Web address. In a 
wavelet image of amplitude or power versus wavelength and 
height, such as we present in Figure 2, a wavetrain of con- 
stant wavelength and amplitude would appear as a constant- 
brightness, vertical band. If the wavelength varied along the 
train, the ridge would tilt or curve. Changes in amplitude 
would appear as varying brightness along the band. 

Following Raynaud et al. (2003, 2004), we used the Morlet 
wavelet with a nondimensional frequency of 6 (wo in Equa- 
tion (1) of Torrence & Compo 1998). The left panel of Fig- 
ure 2 shows the valid region of the temperature profile, for 
which we computed the wavelet amplitudes presented in the 
right panel. Torrence & Compo (1998, their Equation (8)) 
provide a power spectrum "normalization" that allows direct 
comparison to Fourier analyses, which we do below. This ad- 
justs for the wavelet spectrum's geometric spacing between 
frequencies and enables definite integrals of power spectral 
density (PSD) to compute total power in a wavelength range. 
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Figure!. Wavelet amplitudes of the thermal profile. Left: thermal profile from the 1998 November 14 occultation by Saturn (spiky, solid line, bottom axis), 
background profile used for gravity-wave modeling and the "B" case of the global power spectrum, below (smooth, solid line, bottom axis), and Brunt- Vaisala 
frequency derived from background profile (dashed line, top axis). The background profile is a wavelet reconstruction of the thermal profile using only scales 
larger than 60 km. Right: the corresponding unnormalized wavelet amplitude (not power) spectrum, with the same vertical axis. The horizontal axis is the Fourier 
wavelength equivalent to the wavelet scale at each location. The cross-hatched region is the unreliable "cone of influence" (COI) of edge effects that alter the 
amplitudes. The gray level at each point is proportional to the temperature amplitude of a wave at that altitude and wavelength (see scale bar at top, which omits 
levels in the COI for readability). Empirical tests with real and Gaussian noise produced features (at larger periods) at least as strong as the 0.7-K peak at 12.5 
km wavelength, -35 km height. Symbols mark the locations of amplitude peaks fit by gravity-wave models (Table 4). White contours encircle regions of greater 
than 95% confidence per point that the signal stands significantly above the data's global power spectrum. See Section 5 for discussion, including significance 
contours. 



Since amplitude spectra are not spectral densities (i.e., one 
does not integrate them), one must remove the normalization 
to express amplitudes in Kelvins. We thus keep the adjust- 
ment for the power spectra presented herein, but remove it for 
amplitude spectra. This adjustment and the orthogonal basis 
set of the transform allow us to recover accurately the am- 
plitudes of synthetic sinusoidal signals inserted into the input 
data. The cross-hatched region is the so-called cone of in- 
fluence (COI) of the edges of the data. In this region, points 
are close enough to the edge of the data that wavelets at those 
periods extend beyond the data. This effectively averages in 
zeros from outside the data. Structure within the COI is unre- 
liable, so we ignore it. The profile is dominated by structures 
with wavelengths longer than 32 km. Few-kilometer, irregu- 
lar wiggles are superposed, but consistent strength in scales of 
5-30 km is absent. This is reflected in the wavelet transform, 
where amplitudes drop substantially outside the COI. 

We also computed the global wavelet power spectrum (Tor- 
rence & Compo 1998, Equation (22)) by averaging the nor- 
malized wavelet power transform over all valid heights that 
are outside the COI (see Figure 3). This method of computing 
the power spectrum has advantages over both a single Fourier 
transform with a window applied to the data and a similar av- 



erage over a windowed Fourier transform (WET). The WET is 
a local transform similar to a wavelet transform, but its sliding 
window has a fixed width that affects the computation differ- 
ently at each wavelength, introducing wavelength-dependent 
behavior The wavelet transform's window size scales with 
the wavelength and thus has the same effect at each wave- 
length. The local nature of both the WET and wavelet trans- 
form allows one to use the COI to omit points contaminated 
by edge effects. The result is a dramatic reduction in power 
spectrum noise, as shown in Figure 3. The new significance 
test for the global power spectra (see Section 7) could not be 
meaningfully applied without the superior noise rejection of 
wavelets. 

The data's global spectrum shows a downward-curving 
trend with superposed peaks. The general trend appears to 
follow an atmospheric gravity-wave spectrum (see Section 7). 
The peaks could be discrete gravity waves, which we study 
with a model (see Section 6). Finally, we compare the global 
spectra of the data and of our noise tests to establish signifi- 
cance (see Section 7). 

Are the features in the global and local wavelet spectra sig- 
nificant? This is really two separate questions: does the global 
spectrum stand significantly above the noise? and do local 
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Figure 3. Power spectrum of thermal profile (75% cap, B case of normal- 
ized temperature, below) derived in three ways. Thick line: wavelet spec- 
trum averaged outside of COL Medium line: Fourier transform with Hann 
window and 8/3 scaling used, e.g., by Young et al. (2005). Thin line: Fourier 
transform without window or scaling. The wavelet-derived spectrum is much 
smoother than either Fourier-transform-derived spectrum. The non-Hann 
Fourier spectrum rises above the other two curves at short wavelength and 
falls below the Hann version at long wavelength. The high noise level of the 
Fourier-transform-derived spectra would make them unsuitable for use in our 
significance test. 



features stand significantly above the global spectrum? We 
discuss the latter here and present the former in Section 7. 

Torrence & Compo (1998) developed a rigorous signifi- 
cance test based on a lag- 1 autocorrelation noise model, which 
compares the data to the data shifted by one point, on the as- 
sumption of a short correlation length. That test unfortunately 
does not apply to occultation inversions, since inversions are 
correlated from any given depth all the way to the top of the 
atmosphere. Thus, the correlation length varies as much as it 
possibly can, and is short only at the top of the atmosphere. 

We verified that the real and imaginary parts of the power 
transform each had a Gaussian distribution at a given wave- 
length. Then, we followed code comments of Torrence & 
Compo in applying the ^ probability distribution to deter- 
mine the multiplier for the global spectrum that gives the 95% 
confidence level (we assume that if the power is significant at 
a given point, so is the amplitude). In Figure 2, contours encir- 
cle regions with power greater than this level. Note that these 
are not contours of the amplitude transform, since the global 
power spectrum's normalization varies with wavelength. 

The strongest of the significant features that is completely 
outside the COI has vertical wavelength A^. = 12.5 km, a max- 
imum at z = -35 km, and amplitude 0.7 K at that level. This 
feature can be interpreted as a short gravity-wave train, as can 
several others at shorter wavelengths and much lower ampli- 
tudes. It is localized, meaning that it does not extend ver- 
tically over the entire data set. One can see, from the left 
panel in the plot, that this feature lasts at least 1.5 cycles. We 
show in Figure 4 that it extends over four cycles. Young et al. 
(2005) also detected several short wave trains in the Galileo 
Probe data for Jupiter. 
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Aside from their much-lower amplitude at a given wave- 
length, the transforms of the noise data sets are qualitatively 
very similar to that in Figure 2, except that the RS and GS sets 
have decreasing power at lower altitude (as expected). The 
large-amplitude (many-K) oscillations in the noise sets gen- 
erally have wavelengths of at least a scale height, and are in 
the COI. Wavelet amplitude maxima outside the COI for the 
GN set averaged ^0.13 K, and for the three shifts of the RN 
set were 0.19, 0.30, and 0.34 K (these red-noise features ap- 
peared at long wavelengths). Maxima are somewhat smaller 
for the scaled data sets, but this is less meaningful as their 
maxima now come only from the high-altitude region; given 
enough samples we would expect similar maxima. In both 
real and Gaussian cases, about half the tests contained sub- 
stantial, isolated regions outside the COI with amplitudes that 
stood above the surrounding features. These regions were 
broader for the real noise tests, often substantially broader 
than the ones in Figure 2. The maximum transform ampli- 
tude outside the COI often appeared in these features. As we 
note below, most or all of the data transform outside the COI 
stands well above the noise. Again, one assesses this best by 
comparing global spectra. The point here is that one cannot 
tell signal from noise merely by looking at the pattern pre- 
sented. 

Numerous wave features in the literature have amplitudes 
smaller than 1 K, often in noisier data sets than this one (e.g., 
Raynaud et al. 2003, 2004). Such cases likely also have sig- 
nificant, spurious fluctuations in their transforms due to scin- 
tillation noise and other terrestrial atmospheric effects, though 
it is impossible to tell for certain without the test described 
above, as the degree of such noise is rarely reported and time- 
correlated noise could conceivably have been better for those 
observations than that for ours, even if random noise were 
worse. It thus becomes imperative to apply a wavelength- 
dependent significance test, such as that presented above. Had 
we simply used the maximum of our noise transforms, we 
would have rejected all the encircled regions of Figure 2 as 
being below 2a significance. By the same token, it becomes 
difficult to accept such detections without any significance 
test. There are two classes of features, discrete and global. 
The significant, discrete features have contours around them 
in Figure 2. The significant, global features are the bumps in 
the thick line of Figure 3, which we interpret with a gravity- 
wave model in the next section and later compare to the noise 
spectra. 

6. GRAVITY- WAVE MODEL 

The noise transforms do not show the continuous vertical 
ridges of elevated amplitude one would expect of a wave that 
propagates undamped through the entire valid region of the 
temperature profile. However, it may be possible, by fitting 
an appropriate model, to improve sensitivity to such individ- 
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Figure 4. Gravity- wave models. The "Full Profile" panel shows the 75%-cap profile after subtracting the background state (thick line) and the sum of the models 
in the next five panels (thin line). The remaining panels show gravity-wave models fit to these observations and to representative profiles from the (supposedly 
waveless) noise sets. The top of each panel states its source and central reconstruction wavelength (near the wavelength at peak amplitude given in Table 4), with 
"D" representing fits in the real data shown in the leftmost panel. In the model plots, the thick, solid trace shows the reconstructed thermal profile around the 
stated wavelength. The thin, solid trace gives the fitted gravity-wave model. The pair of dashed, horizontal lines shows the altitudes where the data intersect 
the COI; only information between these lines is reliable. Parameters of the "D" wave fits appear in Table 4. Note that, although the RN 20 case has a relatively 
high amplitude for a short stretch in the valid region, its global amplitude (mean amplitude over all valid heights) is much lower than that of the D 20 case. The 
reconstruction for D 20 had a rather narrow range of wavelengths to separate it from D 25, leaving significant power out of D 20. In the left panel, the sum of 
our models does not perfectly fit the data because some models do not fit well at all altitudes and there is a background spectrum of low-intensity activity (see 
Section 7). 



ual gravity-wave modes in the data, using the vertical coher- 
ence that the noise does not exhibit. Raynaud et al. (2003, 
2004) fit such models to identify gravity-wave signatures in 
Jovian temperature profiles derived from two stellar occulta- 
tions. Note that these two papers use different methods. Ray- 
naud et al. (2003) used three chords from one event to de- 
termine both Aj. and the horizontal wavelength. A/,. Through 
their Equation (15) they determine the dissipation level for 
the wave and point out that the observed T peaks are too high 
in the atmosphere to be caused by a gravity wave with the 
derived parameters. In contrast, the event reported by Ray- 
naud et al. (2004) has just a single observation. The horizon- 
tal wavelength is a free parameter and fits are consistent with 
a gravity-wave interpretation. 

We follow Raynaud et al. (2004) very closely to test 
whether waves may be propagating at the wavelengths cor- 
responding to peaks in the global wavelet spectrum shown in 
Figure 3. These correspond to regions outside the COI and at 
Aj < 60 km in Figure 2. The first three columns of Table 4 give 
Az, altitude (Zm-dx), and temperature amplitude {AT[zmdx]) of 
five temperature peaks; these are marked in Figure 2. The re- 
maining columns give the altitude of maximum amplitude of 
the model wave (zo), \, and the wave period. Physically, zo 
represents the altitude where wave damping exceeds the natu- 
ral exponential amplitude growth. The intensity of dissipative 
processes, Aj., and A^ determine this altitude. 

As a wave propagates, variations in the atmosphere's 
steady-state properties affect its amplitude and vertical wave- 
length. In a conservative atmosphere, a vertically propagat- 
ing wave experiences exponential growth of its temperature 
amplitude as a consequence of the exponential decrease of n 
and the requirement for energy conservation. If the amplitude 
becomes large enough for the local temperature gradient to 



exceed the adiabatic lapse rate, F, the wave becomes unsta- 
ble and overturns. Dissipative processes such as molecular 
viscosity, eddy diffusion, thermal conduction, and radiative 
damping can further limit amplitude growth, causing the wave 
to deposit its energy in the background atmosphere (Lindzen 
1981). Furthermore, variations in the background tempera- 
ture and vertical shear in the zonal wind can also cause ob- 
servable changes in the amplitude, vertical wavelength, and 
phase propagation. 

We simulate these effects with a hydrostatic WKB gravity 
wave model in a rotating atmosphere. The model (Matcheva 
& Strobel 1999) includes dissipation by molecular viscosity, 
eddy diffusion, and thermal conduction. We assume that there 
is no vertical gradient in the zonal wind. Calculations that in- 
cluded constant wind shear did not improve the fits. We as- 
sume scales larger than 60 km to represent the atmosphere's 
steady state, although from a single occultation we cannot de- 
termine whether an observed structure is transient or static. 
Our steady-state profile is the inverse wavelet transform of 
scales larger than 60 km in the right panel of Figure 2; it ap- 
pears as the smooth, solid line in the left panel. 

The eddy diffusion coefficient, an input in our wave model, 
parameterizes the intensity of the vertical mixing in the atmo- 
sphere. Published values for Saturn's eddy diffusion coeffi- 
cient differ by more than 2 orders of magnitude (Moses et al. 
2000 and references therein). In our model it varies with al- 
titude and is proportional to n"'/^ (Atreya 1986). We set it to 
500 m"^s at the 0.1-/ibar level to agree with Voyager's UVS 
observations (Smith et al. 1983). The vertical mixing in the at- 
mosphere is important as a wave dissipation mechanism. The 
degree of wave dissipation is determined by the eddy diffu- 
sion coefficient and the wave's vertical and horizontal wave- 
lengths. Since the temperature profile only has vertical in- 
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formation, the horizontal structure is a free parameter The 
exact eddy diffusion coefficient value is thus not critical for 
the wave model because adjusting the horizontal wavelength, 
A/,, within a reasonable range can compensate for reasonable 
changes in the eddy diffusion coefficient without changing the 
vertical structure. 

The wave model parameters are A/,, A^, temperature ampli- 
tude Ar, and phase at the lower boundary. We adjust the 
wave parameters so that the model fits the observations at Zmax 
(see Table 4). Since the background atmosphere's tempera- 
ture varies with altitude, so does \. for each wave. We select a 
short range in A, for each candidate wave, zero all other scales 
in the wavelet transform, and reconstruct the temperature fluc- 
tuations due just to those scales using the inverse transform. 
We compare the simulated wave and the reconstructed tem- 
perature variations within the limits set by the COL The "D" 
panels of Figure 4 present the results. 

For the purpose of wave identification, one would ideally 
like to follow the wave signature over several wavelengths and 
several scale heights, with at least two scale heights below the 
altitude of wave damping. However, the region outside the 
COI becomes progressively shorter for larger wavelengths. 
Therefore, large- A, temperature fluctuations are difficult to in- 
terpret uniquely as signatures of propagating waves. 

We obtain a relatively good fit for the longest wavelength 
that we consider (A^ = 40 km), detecting significant power 
throughout the observed region. The simulated wave has a 
\h = 200 km; the corresponding gravity-wave period is 42 
minutes. The wave exists below the expected altitude of dis- 
sipation (altitude of maximum amplitude), and shows modest 
amplitude growth. Both the amplitude and the phase of the 
simulated wave follow the reconstructed temperature fluctua- 
tions well. However, despite the apparently good fit, we can- 
not rule out alternative interpretations, since the reliable part 
of the reconstruction contains only two wave cycles. 

For the rest of the candidates, we have mixed success in fit- 
ting the observed temperature fluctuations with a single prop- 
agating gravity wave. In general, we are able to fit the recon- 
structed temperature fluctuations (both amplitude and phase) 
well over 2-5 cycles near the peak amplitude, but we are 
not successful in matching the observations throughout the 
entire altitude region. Young et al. (2005) also present wave 
activity in the Galileo Probe data for Jupiter that spans only a 
few cycles before it disappears. As previously (Raynaud et al. 
2004), our requirement for a positive wave identification is 
much more restrictive, since we require a good fit throughout 
the profile's valid region. 

Our wave model fit to the 25 km reconstructed scales does a 
relatively good job following the amplitude of the temperature 
fluctuations throughout the extent of the data, but it fails to fit 
the phase of the observations at the bottom of the occultation 
(z < -70 km). 

In the case of the 20 km structure, we have a good phase 
match throughout the sampled region, but the amplitude ex- 
hibits a double peak that is difficult to explain with a simple 
gravity-wave model, with or without a constant vertical wind 
shear. The amplitude fits well within four wavelengths of the 
top and bottom of the valid region, but overestimates the ob- 
served fluctuations in the central part of the profile. A variable 
wind shear might explain the observed amplitude. 

The 12.5 km structure has a well defined power maximum 
that is nicely separated from the COI. This and the relatively 
small vertical scale make the wave analysis more robust. We 
achieve a good amplitude and phase fit for four wavelengths 



for z > -50 km. At z = -50 km, the amplitude shows a 
very fast increase with altitude, which is inconsistent with our 
model. Such a fast amplitude variation can be an indication 
of significant vertical shear in the zonal wind. If this were the 
case, however, we should also detect a sudden change in A,., 
but the observed change is not very large. 

The temperature reconstruction of the 6.5 km candidate 
shows a significant variation in amplitude. These variations 
can be modeled with some success if one assumes that the 
background temperature profile retains scales, L^, below 60 
km (e.g., L, > 30 km). For our background temperature pro- 
file (L, > 60 km), the model fits well in both amplitude and 
phase for 4.5 wave cycles at altitudes below -100 km. How- 
ever, the model does not reproduce the multiple amplitude 
peaks. The significant amplitude variation suggests beating 
wave modes whose wavelengths are unresolved by the data. 
The observed amplitude can be well modeled by superposing 
two sine waves that have vertical wavelengths 5.96 km and 
6.66 km respectively. 

For a given eddy diffusion coefficient, K, waves with short 
Xy and long A/, dissipate low in the atmosphere. This places 
a natural filtering mechanism for wavelengths that can prop- 
agate at a given altitude. In order for a wave with a short 
Aj. (e.g., 6 km) to propagate at the probed pressure levels. A/, 
must be rather small (^10 km). The wave then comes close 
to violating the hydrostatic approximation (A^. << A/,). It also 
raises a question about the detectability of such a wave by a 
stellar occultation, which averages the properties of the atmo- 
sphere along the line of sig ht. The length of this averaging 
on Saturn is about VlnrH = 4500 km, where r is the plan- 
etary radius. Sicardy et al. (1999) demonstrate that the am- 
plitude of light-curve fluctuations resulting from the presence 
of a monochromatic wave with a projection of the horizontal 
wavelength along the line of sight / = A/,/cos0; is not signifi- 
cantly reduced if 

where 0° < 6*/ < 90° is the angle between the line of sight and 
the horizontal wave propagation direction. For Saturn, 77 = 18. 
If the 6.5 km feature is to be interpreted as a gravity wave 
with a horizontal wavelength A/, = 10 km, the wave must be 
propagating at an angle 61 > 85° in order to be detected dur- 
ing the occultation. In other words, we would have detected 
a wave that propagated almost along the planetary meridian. 
In this respect, of all the waves that might be present in the 
atmosphere, we preferentially detect waves that propagate in 
the horizontal at large angles 61. In summary, the arguments 
for a wave interpretation of the very short scales (5-8 km) 
present in the temperature profile are easiest to accept if the 
wave field is isotropic. 

To see how well the model-fitting method rejects noise, we 
performed the same analysis on one randomly selected profile 
from each noise set; results are in the final eight panels of Fig- 
ure 4. For each type of noise realization we reconstruct two 
scales, one short (less than 30 km) and one long (about 40 - 
60 km), and fit gravity-wave models. The fits to reconstruc- 
tions of noise and data show similar qualities. Short scales fit 
well over a few wavelengths but not through the entire alti- 
tude range, and large scales fit relatively well throughout the 
vertical range, but the region outside the COI contains less 
than two full wave cycles. This makes it difficult to identify 
an atmospheric wave based only on the goodness of fit. How- 
ever, there are ways to distinguish some waves from noise. 
First, our detected waves have much higher amplitude than 
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Figures. Normalized temperature profiles derived from the 75%-cap in- 
version. The thick profile uses the wavelet-smoothed background profile of 
Figure 2 (B case). The thin profile uses the mean of the temperature profile 
in the altitude range presented (M case). 

the noise at the same scale. The relevant scales appear signif- 
icant in the global wavelet power spectrum (Figure 3), even 
though most span only a short vertical range. Second, in the 
noise fits most of the power (including the peak) is typically 
within the COl, and the amplitude varies with altitude faster 
than expected for a wave. Unfortunately, the latter are tenden- 
cies rather than robust discriminators; real waves can do the 
same, making it difficult to specify robust criteria to discrimi- 
nate between waves and noise-induced features. For example 
the fits of the 40-km data structure (D40) and the 40-km noise 
structure (GN40) look very similar. The amplitude of the fit 
to the real data is, however, more than five times larger than 
the fit to the noise. This signal-to-noise assessment is thus the 
only reliable discriminator known to us. 

We conclude that, while this method may do well at quanti- 
fying the properties of strong waves, it is not a good discrim- 
inator between weak waves and noise. A number of criteria 
should be satisfied for the entire region outside the COl to 
make a positive wave identification: 

1 . the wave amplitude should be many times the mean am- 
plitude of the worst-case real noise at that A^; 

2. a wave model should fit both the amplitude and phase 
of the observed temperature fluctuations well; 

3. the structure should have more than a few cycles. 

Note that structures failing one or more criteria may still be 
(or contain) waves. We discuss the broad spectrum of weaker 
gravity waves in the next section. 

7. POWER SPECTRUM 

To compare the power spectrum of temperature fluctua- 
tions in Saturn's atmosphere with that in other atmospheres, 
we calculate global_ spectra for normalized thermal profiles, 
(T - f)/f, where f is the mean temperature. While others 
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Figure 6. Power spectra of inverted noise-test light curves (thick lines), and 
power-law fits (thin lines). The words "max" and "mean" indicate global 
spectra generated by those two methods from each of the four types of noise 
tested (see the text). Parameters of these fits appear in Table 3. We compare 
them to the data in Figure 7. 



(e.g.. Young et al. 2005; Allen & Vincent 1995) have sim- 
ply chosen regions of their data where the background was 
clearly isothermal to calculate f , we do not have this luxury, 
so we computed it in two ways. In the first, f is the mean 
temperature j^l 48. 04 K) in the good region of the profile. In 
the second, f is the wavelet-reconstructed thermal profile for 
\ > 60 km (the smooth, solid line in the left panel of Figure 
2). The choice of 60 km is arbitrary, hence the two cases. We 
call them M and B, for mean and background, respectively, 
and present them in Figure 5. 

From their wavelet transforms, we derived the global power 
spectrum for both normalized temperatures from our nominal 
profile (75% cap), as described in Section 5. To derive noise 
spectra, we calculated the wavelet spectrum of each noise in- 
version, then found the mean power at each wavelength over 
all the altitudes and inversions in each set, excluding points 
in the COL Since we demonstrated above (Figure 4) that real 
noise can induce wave-like features, we also derived spectra 
based on the maximum rather than the mean, to ensure that 
anything we call real must stand well above these spurious 
features. Figures 6 and 7 and Table 3 present the resulting 
power spectral densities and power-law fits, computed with 
the full, wavelet-normalized power spectrum for comparison 
to other work. The table also gives power-law fits to tempera- 
ture amplitude, with neither wavelet nor temperature normal- 
izations applied. 

The noise sets' power spectra follow an m~^ dependence, 
where m is the wavenumber, that is very consistent along their 
entire lengths, save for small tails near the Nyquist frequency. 
The tails are well below the resolution of the uninterpolated 
light curve, so we exclude them from the fits. The non-white 
noise should be due almost entirely to Earth's atmosphere. 
The m~^ power law is unfortunate, as several atmospheric pro- 
cesses share it. These include the short-scale end of the uni- 
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Figure 7. Global power spectrum of the data (thick line) compared to data fits (smoothly curving lines, see key) and noise spectrum fits (straight lines). This 
compact presentation relates the main points of this paper; see Section 7 for interpretation. Mean-normalized data and fits are on the left, background-normalized 
on the right. The noise spectra fits are error estimates for the data spectrum under our eight noise models for power spectra of inversions. Table 3 gives parameters 
of these fits. The three characters of each noise label indicate, in order, the type of noise set (Real or Gaussian), whether it is a Normal or Scaled set, and whether 
the spectrum came from a vertical average over the maXima or Means of the power at each location in the wavelet spectra of the inversions in that noise set. For 
example, "RNX" means Real noise. Normal (non-scaled), averaged vertically over the maXima. Mean spectra have dashed lines; maximum spectra have dotted 
Unes. Standing many times higher than the linear noise spectra, the data show peaks superposed on the modified Desaubies function (Equation (2)) that describes 
the "universal spectrum" of gravity waves for Earth. The four smoothly curving traces are the modified Desaubies fits given in the key and Table 5. The peaks 
are the discrete waves shown in Figure 4, analyzed in Section 6, and presented in Table 4. They cause the large reduced-x^ values of Table 5. 



versal gravity-wave spectrum observed on Earth and Jupiter 
(see below) and short-scale, two-dimensional (2D) turbulence 
(Harrington et al. 1996; Travis 1978), although 2D dynam- 
ics, at least, likely break down at scales larger than are rel- 
evant to the current observations. The scaled noise sets had 
lower power (see Figure 6 and constants in Table 3), but oth- 
erwise behaved similarly to their unsealed brethren. The spec- 
tra computed by taking the maximum lay a factor of 4 or more 
above those computed by taking the mean, which the fit con- 
stants of Table 3 reflect. 

Figure 7 presents the data and fits to the noise spectra (and 
spectrum models; see below). At all wavelengths, the power 
spectrum of the data is 2 - 10 times higher than even the 
most conservative noise spectrum (labeled "RNX") and de- 
parts from the very linear shape of the noise spectra. We can 
confidently accept that the data's spectrum reflects properties 
of Saturn's atmosphere rather than noise. The log plot allows 
by-eye S/N calculation under different noise models. For ex- 
ample, at a wavelength of 32 km, the data PSD is lO"' and the 
RNX noise spectrum is 2x lO"'*, indicating a 5cr detection of 
the planetary spectrum there. Under GNX noise it would be 
lOcr, and under GNM noise it would be a SOcr detection. We 
thus conclude that comparison to real, not Gaussian, noise 
is necessary and that averaging over the maximum power at 
each wavelength and location in the noise set (rather than the 
mean) is necessary because of the spurious, wave-like features 
induced by real noise. 

We can now investigate the spectrum of gravity waves on 
Saturn, though with the caveat that temperature profiles de- 
rived from occultations give an average over a long horizontal 
path through the atmosphere (see above). On Earth, the spec- 



trum of temperature fluctuations is constant enough over sea- 
son and location to be called the "universal spectrum" (Van 
Zandt 1982; Balsley & Carter 1982; Dewan et al. 1984; Vin- 
cent 1984; Smith et al. 1987; Tsuda et al. 1989). It has been 
observed on Jupiter (Young et al. 2005) as well, and follows 
the modified Desaubies function. 



p{m) = a^r— 



(m/m^,y 



(2) 



g^m^ 1 + {m/m^y^' 
where p is the PSD, a is a unitless constant, A^ is the 
Brunt- Vaisala frequency, m* = |r|/2(Tj' = 27r/L, is the criti- 
cal wavenumber, L* is the critical wavelength, aj is the rms 
temperature fluctuation of {T - f), s is the long-wavelength 
power-law exponent, and t is the short-wavelength exponent. 

The characteristic wavenumber m, divides the spectrum 
into two regions: a small wavenumber region (m < m,) and a 
large wavenumber region (m > m*). The shape and magnitude 
in the small-wavenumber region (large vertical wavelengths) 
is believed to be dominated by the wave source characteristics 
and is proportional to m\ The exponent s for the terrestrial at- 
mosphere is not well constrained, though it is typically about 
1 . In the large-wavenumber regime (small wavelengths) satu- 
ration and/or dissipation processes are believed to control the 
wave spectrum (Gardner 1996 and references therein). In this 
region the spectrum is proportional to m~', where f = 3. The 
overall amplitude of the spectrum is controlled by the dimen- 
sionless parameter a, which may represent the wave genera- 
tion mechanism. 

Young et al. (2005) applied terrestrial parameters a = 0.1, 
s = 0, and f = 3, and derived A^ = 0.0176 s ^ F = 2.1 1 K km', 
and (77 = 5.0 K to compute L* = 30.3 km and thus the De- 
saubies function for Jupiter. Note that in their paper, a typo- 
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Table 5 

Parameters of the Modified Desaubies Fits 



Case 


a 


(km"l) 


s 


t 


(km) 


X^/dof*' 


M Jupiter 


0.1000 


0.207 


0.00 


3.00 


30.3 


780 


B Jupiter 


0.1000 


0.207 


0.00 


3.00 


30.3 


791 


M nominal 


0.1000 


0.273 


0.00 


3.00 


23.0 


720 


B nominal 


0.1000 


0.367 


0.00 


3.00 


17.1 


681 


Mfit4 


0.1098 ±0.0016 


0.564 ± 0.043 


-1.71 ± 0.09 


3.92 ± 0.06 


11.1 ±0.9 


17 


B fit4 


0.1089 ±0.0015 


0.725 ± 0.056 


-1.95 ± 0.07 


4.23 ± 0.08 


8.7 ± 0.7 


13 


Mfit3 


0.0848 ± 0.0015 


0.299 ± 0.006 


0.00 


3.49 ± 0.02 


21.0 ±0.4 


21 


B fits 


0.0872 ± 0.0016 


0.313 ±0.006 


0.00 


3.58 ± 0.02 


20.1 ±0.4 


19 



^L*=27r/m,. 

*" dof = degrees of freedom. Note that we fit only the broadband spectrum, not the waves. 



graphical error gives m, = Taj /2, but the value is calculated 
correctly. Their Equation (3) gives a curve that lies a factor of 
^ 2tt lower than the curve presented in their Figure 7, using 
the parameters given above, though the shape is the same. It 
is unclear whether the plotted data have been similarly shifted 
or whether the two traces do not in fact follow one another. 

Figure 7 presents Equation (2) evaluated for the parameters 
given in Table 5. Some of these curves have free parameters 
fit with a Levenberg-Marquardt minimizer; the uncertainties 
on those parameters are from the minimizer's covariance ma- 
trix. The "Jupiter" case in Table 5 and Figure 7 uses m, from 
Young et al. (2005), but A^ = 0.01 165 s' and T = 0.992 Kkm', 
appropriate to our data. The poor visual fit and high x^ show 
that it misses the data by a substantial margin. Using crj- calcu- 
lated from the B and M cases (1.18 and 1.98 K, respectively) 
yields the m, and L* given in the "nominal" cases (i.e., Satur- 
nian values and no free parameters). Things are better, but not 
good. However, the long atmospheric path of the occultation 
ray likely reduced aj well below what would be observed in- 
situ; for the Galileo Probe data, ctt- = 5.0 K, for example. For 
the case where all parameters vary, a is still comfortably close 
to 0.1. Fixing s = puts a closer to 0.085. In all cases with 
free exponents, the exponents move substantially away from 
their nominal values. In particular, f is now decidedly steeper 
than the -3 power law usually cited for gravity wave spectra. 
Since waves with small A, are more likely to have small A/,, 
the steeper power law may be due to preferential filtering of 
short wavelengths due to averaging along the occultation ray 
path. 

The reduced x^ is still high, indicating that the undulations 
in the spectrum are significant: there is more here than a back- 
ground spectrum of gravity waves. The gravity- wave model 
presented above fits averages over the sections of the wavelet 
transform that produce these bumps. That model assumes a 
single wavelength per wave, so we cannot use its output to 
make a complete, simultaneous model of background and dis- 
crete waves that might fit the global spectrum more perfectly. 

8. DISCUSSION OF THE SPECTRUM 

The one-dimensional power spectrum of terrestrial grav- 
ity waves exhibits a nearly universal behavior in its large 
wavenumber region (see the previous section). This state- 
ment is based on a large number of observations of the verti- 
cal structure of the horizontal wind, temperature, and density 
fluctuations of the present wave modes using diverse experi- 
mental techniques. The shape of the spectrum at large verti- 
cal wavenumbers (short vertical wavelengths) is consistently 
independent of time, place, and altitude. This is usually at- 



tributed to a saturation process that limits the wave amplitude 
growth. 

Theoretical work diverges on the nature of the satura- 
tion process at work. The existing theories use different 
physical mechanisms for dissipating wave energy including 
shear and convective instabilities (Dewan & Good 1986), cas- 
cade processes (Dewan 1991), wave-induced Doppler effects 
(Hines 1991), and wave-induced diffusion (Gardner 1994; 
Zhu 1994). The different theories all predict the same shape 
and behavior for the PSD, including the -3 slope of the 
large wavenumber tail and the presence of a characteristic 
wavenumber m* that limits the range of the saturated spec- 
trum. Since m, oc 1 /aj, the theories predict a decrease in m* 
with altitude as m* ex e"'/'*'^, and the terrestrial data agree. 

Despite the diversity of wave-generating mechanisms and 
atmospheric thermal structure on different planets, the satura- 
tion theories are based on rather general concepts, implying 
similar spectra on different planets. Confirming this on other 
planets is challenging as it requires local measurements at a 
variety of altitudes and horizontal locations on each planet, 
yet there have been only a few atmospheric entry probes. Oc- 
cultations average over a long path length, underestimating 
the wave amplitudes. This effect depends on wavelength and 
probably modifies the true spectral slope, especially at large 
wavenumbers. Indeed, in Earth's atmosphere a typical value 
for / is 2.5 - 3 (Allen & Vincent 1995), whereas our best fits 
show a slope of 3.5 - 4.2. The higher values for t might also be 
real, reflecting actual differences between the dominant satu- 
ration processes acting in Earth's and Saturn's atmospheres. 
Gravity wave studies in Earth's atmosphere are based on thou- 
sands of temperature profiles taken at different locations, al- 
titudes, and times of the year, using ground-based, airborne, 
and satellite-based techniques, whereas we are discussing a 
single temperature profile. On the other hand, the value of 
L* ss 15-25 km, based on the observed temperature variance 
(see Equation (2)), is in good agreement with the value de- 
termined using the Desaubies fits (see three-parameter fits in 
Table 5). 

Where should we focus our efforts for wave detection? In 
addition to seeking individual wave modes that may reveal 
information about atmospheric structure before losing their 
identity through interactions, one can look at the spectral char- 
acteristics at small wavenumbers (large vertical wavelengths). 
The theory says that at wavenumbers smaller than m,, the 
waves are not affected by the saturation processes (whatever 
they are) and retain the spectral characteristics of the gen- 
eration mechanism. Numerical simulations of convectively 
generated gravity waves show that there is also a dependence 
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between the characteristic vertical wavenumber m, and the 
depth of the convective cell generating the waves (Fritts & 
Alexander 2003). 

9. CONCLUSIONS 

We have analyzed a light curve based on IRTF observations 
of the 1998 November 14 occultation of GSC 0622-00345 by 
Saturn (Paper 1). We presented and analyzed isothermal light- 
curve fits, an atmospheric thermal profile, a wavelet analysis, 
gravity wave modeling, and power spectra. The derived ther- 
mal profile varies over 142 - 151 K in the pressure range 1 - 
60 /ibar at a latitude of 55.°5 S. The vertical temperature gra- 
dient is removed by more than 0.2 K km ' from the adiabatic 
lapse rate, indicating that the stratospheric region sounded by 
the occultation is statically stable. Our thermal gradient pro- 
file shows the same alternating-rounded-spiked appearance of 
other occultation profiles, including one for Saturn (Cooray 
et al. 1998). This shape has previously been interpreted as ev- 
idence of gravity wave breaking (Raynaud et al. 2003, 2004; 
Young et al. 2005). 

Our new noise tests, based on real noise sampled from our 
light curve's upper baseline rather than synthetic, uncorrected 
Gaussian noise, showed that atmospheric scintillation (and 
similar correlated noise sources such as spacecraft pointing 
drifts) can introduce relatively strong, spurious temperature 
fluctuations into the thermal profile derived by inverting an 
atmospheric occultation light curve. For our data, the effect 
was many K for the raw inversions, but it mainly appeared 
at longer wavelengths in power spectra. However, the ampli- 
tudes at the shorter wavelengths found outside the COl were 
still of order 0. 1 - 1 K, comparable to the amplitudes of the 
gravity waves often identified in occultation data sets. We thus 
developed several significance tests for our power-spectrum 
analyses, and note that without such tests, one must be skep- 
tical of gravity-wave detection claims in ground-based (and 
possibly some space-based) occultation inversions. 

We used a wavelet analysis to search for localized gravity 
wave trains. Based on the wavelet power and the use of sig- 
nificance tests only, the strongest candidate had an amplitude 
of 0.7 K, Aj. = 12.5 km at Zmax = -35km, and lasted four cy- 
cles. It and several shorter-wavelength features stand over the 
global power spectrum at the 95% confidence level or greater. 
Alternative explanations for the observed periodic structures 
include sound waves, planetary waves, and non-transient fea- 
tures. Without knowing the horizontal structure of the wave 
we cannot rule out sound waves or planetary waves as the 
cause for the temperature fluctuations. However, one can 
make the argument that, for a given A,, planetary waves typ- 
ically have lower (time) frequencies and therefore dissipate 
lower in the atmosphere than the valid region of our profile. 

To take into account both phase and amplitude informa- 
tion, which could improve sensitivity over our amplitude- 
based noise limit, we fit a gravity-wave model to the ampli- 
tude and phase of the strongest features in the valid region of 
the wavelet spectrum, following Raynaud et al. (2004). The 
model calculates the temperature amplitude vs. height of a 
single, damped wave mode propagating throughout our pro- 
file's valid region. We performed fits both with and without 
a new parameter for constant vertical wind shear. The added 
parameter did not improve the fits, so we report the shearless 
fits here for consistency with prior results from this model. 

Our best candidate for a gravity wave that propagates con- 
tinuously through the valid region, as assessed by our model, 
had Aj; = 40 km, Xh = 200 km, and a period of 42 minutes. 



At this altitude the wave is not strongly affected by dissipa- 
tion and achieves a maximum amplitude above the observed 
atmospheric region exceeding 1 K. According to Equation 
(1), this wave was detectable if the angle between the line 
of sight and the horizontal direction of wave propagation ex- 
ceeded 74°. However, similar fits to quasi-periodic features in 
thermal profiles derived from isothermal (i.e., waveless) light 
curves with real noise gave some fits of similar appearance. 
This demonstrates the need for a study of real noise to estab- 
lish an amplitude criterion that discriminates real waves from 
noise. The reason wavelike (i.e., sinusoidal) features arise 
out of the noise is simple: by reconstructing only a limited 
range of wavelengths, the resulting profile is certainly sinu- 
soidal and has a favored period. Wave models are sinusoidal, 
and ours even has parameters that allow the phase and am- 
plitude to vary, so we will get a good fit if the amplitude of 
the reconstructed data does not vary much. We find that this 
circumstance occurs in the noise data sets often enough to re- 
quire at least criteria for significant wave amplitude. Criteria 
involving the number of cycles outside the COl or the phase of 
the reconstructed data could potentially provide even-more- 
stringent limits. 

We derived global power spectra from our wavelet trans- 
forms, using only data outside the COl. This method greatly 
reduced the noise level of the spectra, which stand everywhere 
2-10 times above the noise level calculated using real-noise- 
contaminated isothermal light curves. The power spectra fol- 
low the modified-Desaubies form of the universal spectrum 
of gravity waves, though with a slightly more negative high- 
wavenumber exponent. Superposed on this spectrum one sees 
the signature of the discrete wave structures discussed above. 
That we see both the universal spectrum and individual fea- 
tures fit well by gravity-wave models lends confidence that we 
are indeed looking at a signal dominated by gravity waves. 

The amplitudes of all five wave-like features that we ana- 
lyzed are well above the noise level. The wave model used 
to fit these features is based on the assumption that the waves 
propagate independently of each other Interactions between 
the wave modes might be able to explain the discrepancies 
between the observed wave fluctuations and the model's sin- 
gle gravity waves. Note that the discrepancies are more sig- 
nificant when derived wavenumbers are larger than the ob- 
servationally derived characteristic wavenumber. This is the 
saturated part of the spectrum where wave-wave interactions 
determine wave behavior and vertical propagation. 

Direct comparison of our derived temperatures with pre- 
vious Saturn occultation measurements (see, e.g.. Table VII 
of Hubbard et al. 1997) requires care, both because of dif- 
ferences in the assumed mean molecular mass and because 
Saturn's mean stratospheric temperature is strongly affected 
by seasonally varying insolation (Bezard & Gautier 1985). 
Heating by inertia-gravity waves might also be important, 
at least in some regions and/or seasons (Cooray et al. 1998; 
French & Gierasch 1974; Young et al. 1997). Detailed mod- 
eling of stratospheric temperatures requires taking account of 
non-LTE effects as well (Appleby 1990). A clearer picture of 
zonal and seasonal variations in Saturn's stratospheric struc- 
ture should emerge when Earth-based stellar occultations can 
be viewed in the context of data from the Cassini orbiter. 
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